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ABSTRACT 

We present an analysis of the variations seen in the dispersion measures (DMs) of 20 
millisecond pulsars observed as part of the Parkes Pulsar Timing Array project. We 
carry out a statistically rigorous structure function analysis for each pulsar and show 
that the variations seen for most pulsars are consistent with those expected for an in- 
terstellar medium characterised by a Kolmogorov turbulence spectrum. The structure 
functions for PSRs J1045-4509 and J1909-3744 provide the first clear evidence for 
a large inner scale, possibly due to ion-neutral damping. We also show the effect of 
the solar wind on the DMs and show that the simple models presently implemented 
into pulsar timing packages cannot reliably correct for this effect. For the first time 
we clearly show how DM variations affect pulsar timing residuals and how they can 
be corrected in order to obtain the highest possible timing precision. Even with our 
presently limited data span, the residuals (and all parameters derived from the tim- 
ing) for six of our pulsars have been significantly improved by correcting for the DM 
variations. 
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1 INTRODUCTION 

The Parkes Pulsar Timing Array (PPTA) is a project which 
aims to take advantage of the extraordinary rotational sta- 
bility of short period (millisecond) radio pulsars. The prin- 
cipal goal of the PPTA is to make a direct detecti on of grav- 
itational waves (|Hobbsl 120051 : iManchesterl [2006'). For this 
purpose it is necessary to measure weekly times of arrival 
(TQAs ) of ~20 pulsars with a precision between 100 and 
500 ns (| Jenet et al.ll2005h . In order to achieve this goal all 
systematic errors in the TOAs must be considered and, if 
possible, corrected. One such correction is the delay caused 
by the plasma between the pulsar and the Earth. Most of 
this plasma contribution is from the interstellar medium, 
but the contribution of the solar wind cannot be neglected 
and the ionosphere will occasionally be important. The dis- 
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persion in the plasma is a linear effect and can, in principle 
be corrected exactly. The group delay, tg(u), is related to 
the integral of the electron density, n e , from the Earth to 
the pulsar, t g (u) — DM/(AV 2 ), where 

DM = f n e dl, (1) 
Jo 

is the "dispersion measure". The dispersion constant K = 
2.410 x 10~ 4 MHz" 2 cm' 3 pes -1 , v is the observing fre- 
quency and L the distance from the Earth to the pulsar. 
When TOAs, t g i and t g 2, are measured at two frequencies, 
v\ and V2, the DM can be estimated using 

DM = K ta _l - jgL . (2) 
u 2 - v x 

A rough estimate of the DM of a (radio) pulsar is gener- 
ally obtained from the discovery observations. This estimate 
can be quickly refined by re-observing the pulsar with more 
widely separated frequencies. Measured DMs for currently 
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known radio pulsars lie between 2.38 and 1235 cm 3 pc and 
from 2.65 to 244cm~ 3 pc for the subset of millisecond pul- 
sars (jManchester et al.ll2005l fl. 

Precise measurements of DM show that it often has sig- 
nificant time variations. A time delay of 100 ns at an observ- 
ing frequency of 1400 MHz, the accuracy goal of the PPTA, 
corresponds to a DM variation of 4.72 x 10 -5 cm -3 pc. Vari- 
ations of this order can occur in the ionosphere only for 
zenith angles in excess of 80° or during major geomagnetic 
storms, so ionospheric corrections will not normally be nec- 
essary. At this level of timing precision, significant variations 
in DM can occur due to the solar wind even when the pul- 
sar is 60° away from the Sun. Variations in the interstellar 
plasma DM result from plasma turbulence and usually have 
a Kolmogorov power spectrum, implying that the variations 
are larger over longer time-scales. In the pulsars observed 
by the PPTA project, such DM fluctuations can reach levels 
that require correction within a few days or weeks. 

It is clear that the goals of the PPTA project cannot 
be reached without measuring the DM and correcting for 
the plasma delay for each observation. The most precise 
TOA measurements are usually obtained at a frequency of 
1400 MHz, but the only dual-band receiver available at the 
Parkes telescope is at 685 and 3100 MHz. Thus the DM vari- 
ations are measured using the dual-band system at different 
times than the TOA observations at 1400 MHz. Since the 
DM varies relatively smoothly, the DM correction can be 
interpolated to the epoch of the primary TOA observation. 
In this paper we use the first few years of DM measurements 
to test methods of correcting for the solar wind, to study the 
interstellar plasma turbulence and to derive algorithms for 
correcting the TOA measurements. 



2 CAUSES OF DM VARIATIONS 

The contributions of the ionosphere and the solar wind have 
been well-studied and can be estimated by various methods 
independently of the PPTA. The total electron content of 
the ionosphere ("TEC") is regularly monitored because it 
is needed to correct the Global Positioning System (GPS) 
navigational system. A monitor is located at the Parkes ob- 
servatory, but it is seldom necessary to make this correction. 
Corrections for the solar wind are implemente d in the stan- 
dard pulsar timing codes t empo and tempo2 (|Hobbs et al.l 
2006; Edwards et al.l 2006). However, these assume a spher- 
ically symmetric solar wind with a constant scale factor and 
do not model observed variations in wind density with lat- 
itude, longitude and time which can be as much as a fac- 
tor of four at any distance (|McComas et alj|2000h . The ear- 
lier package, tempo, assumes a higher density compared to 
tempo2. Neither of these is adequate for the desired PPTA 
precision. 

The DM due to the interstellar medium varies for a 
variety of reasons. For example, variations are known to 
occur for some pulsars within supernova remnants, when 
wisps of ionised gas drift across the line of sight to the pul- 
sar. For instance, the DM of the Vela pulsa r decreased at a 
rate of 0.04 cm -3 pc yr" 1 from 1970 to 1985 ^Hamilton et all 
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Il985h . Similarly the Crab p ulsar shows varia tions up to 
0.02cm" 3 pcyr _1 over 15yr (|Lvne et all Il988f ). Pulsars in 
binary systems which exhibit eclipses show DM variations 
from the ionised envelope of the companion object. These 
have been measured for t wo of the binary p ulsars in the glob- 
ular cluster 47 Tucanae (|Freire et alj|2003f ). The DM change 
of 0.0065 cm" 3 pc for one of these, PSR J0023-7203J, is 
100 times the level that would require correction for the 
PPTA pulsars. Even larger changes have been observed in 
PSR B1259— 63 which is in orbit with a massive B2e star, 
reaching 10.7 and 7.7 cm -3 pc d uring the periastro n passages 
of 1994 and 1997, respectively (|Wang et alj|2004h . 

The DM also varies due to turbulent spatial variations 
which drift across the line of sight between the Earth and the 
pulsar. These have commonly been characterised in the liter- 
ature as linear slopes in DM. Measurem ents of such dDM/d t 
values for four pulsars were discussed bv lBacker et al.l (|l993l ) 
who proposed that dDM/dt oc (DM) 1/2 and modelled the 
variation using wedges of enh anced density. Obs ervations of 
374 pulsars were presented bv lHobbs et al, (|2004l ) who found 
the same relationship. However, a better characterisation of 
the DM variations can be made using the theoretical spa- 
tial characteristics of a turbulent process. As shown later, 
in a turbulent model the relation dDM/dt oc (DM) 1 / 2 arises 
naturally and does not require a wedge model. Th e spatial 
power spectrum of electron density was defined by iRickettl 
\ 19901 ) as 



P(q) = C 2 n q- 



2-k/1 <q< 2n/U, 



(3) 



where C 2 scales the power spectrum (and thus the total en- 
ergy in the process), /3 is the power- law exponent (which is 
11/3 for a Kolmogorov spectrum), l is the outer scale and 
li is the inner scale. Physically the outer scale is identified 
with the largest scale in the medium, typically the size at 
which it becomes inhomogeneous, and the inner scale is the 
scale at which dissipation occurs. Energy is introduced at 
some scale between l and U, supporting the spectrum. This 
energy 'cascades' in frequency to h where it is dissipated. 
Turbulent variations in electron density can be estimated 
from DM variations and various diffractive phenomena such 
as angular scattering, pulse broadening and intensity scin- 
tillations. Diffractive variations are caused by much smaller 
scale fluctuations in density and thus probe different regions 
of the spatial spectrum than DM variations. Diffractive vari- 
ations are modulated by refractive variations which can be 
used to probe scales between the diffractive and the DM 
scales. All these observed variations result from the motion 
of the line of sight through the scattering medium, thus spa- 
tial variations of scale s are associated with temporal varia- 
tions of scale T by s = VT where V is the velocity of the line 
of sight with respect to the scattering plasma. Therefore, the 
inner time-scale ti corresponds to U = Vn. 

Radio scattering observations, such as those mentioned 
above, are directly sensitive to a statistic called the "phase 
structure function", D${t), which is defined by 



D+( T ) = {[4>(t + T)-<Kt)] a ), 



(4) 



where (j> is the geometrical phase delay between the source 
and the observer and the angle brackets denote an ensemble 
average. For a power-law spatial spectrum with 2 < j3 < 4 
between the inner and outer scales, the structure function 
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L> is given by (|Rickettlll977h : 

D,(r) = (r/r d r , 



(•») 



where a = /3— 2. Similarly, we can define a structure function 
for the DM variations: 



Bdm(t) = ([DM(t + t) - DM(f)f). 



(6) 



At scales that are larger than the scale of refractive scintil- 
lation, the geome trical phase approaches the actual phase 
l|Coles et al.lll99l1 ) and these two structure functions can be 
related through the dispersion relation (Equation [2]| , 

Bdm(t) = (Ku/2w) 2 D^t). (7) 

The structure function was first used to investigate 
DM variations bv iRickettl (| 198ST l . Subsequently, the tech- 
niqu e was applie d to P SR J1939+213 4 (B19 37+21) 
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Kaspi ct al 
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iRamachandran et alJ 
that the structure function was consistent with a power-law 
fluctuation spectrum wi th index (5 betw een 11/3, the 
Kolmogorov value, and 4. iKaspi et al.l (jl994) continued this 
work and obtained f3 — 3.874±0.011. From the approximate 
agreement of the diffractive timescale Td co mputed from 
Equat ion [5] and the directly measured value, ICordes et al.l 
(1990) inferred that the inner scale of the fluctuation 
spectrum, U, was le s s than about 2 x 10 7 m. Recently, 
IRamachandran et al.l (|2006l ) extended the data-span to 
20 yr and obtained /3 = 3.66 ± 0.04 which is consistent 
with the value expected f or a Kolmogorov spectrum and 
suggested ~ 1.3 x 10 9 m. ICognard fc Lestradei (|l997T l pre- 
sented the DM variations of a different millisecond pulsar, 
PSR J1824-2452 (B1821-24), and obtained = 3.7 ± 0.2 
which is also consistent with a Kolmogorov spectrum. Dis- 
persion measure varia t ions w ere measured for six pulsars by 
IPhillips fc Wolszczanl (Il99ll ) and structure functions were 
obtained for PSRs B0834+06, B0823+26 and B0919+06. 
Measured power-law indices were in the range 3.77 to 3.87 
with uncertainties of 0.04 or less. 

Assuming that the DM variations are due to turbu- 
lence then, from the definition of the structure function, the 
"slope" dDM/dt, averaged over an interval t, will be a ran- 
dom variable with an rms of [ J D DM (r)] 1/2 /r. This can be 
related to the mean DM value by the distance to the pulsar 
L as both DM and Ddm(t) are proportional to L. Thus the 
observed result that dDM/dt oc (DM) 1 ' 2 is expected for any 
turbulent medium an d does not require a d hoc models such 
as the wedge model of lBacker et alJ (| 1993T ) . This proportion- 
ality will be valid for spatial scales Vt that are less than the 
parsec scale of interstellar clouds, since it assumes that con- 
tributions to the DM fluctuations from various points on the 
line of sight add incoherently. 

There have been three dissipation mechanisms dis- 
cussed in the literature: ion cyclotron damping (which is the 
primary mechanism in the solar wind) , Landau damping and 
ion-neutral collisional damping. It is not thought that Lan- 
dau damping is important in the ISM (jMinter fc Spanglerl 
1 19971 ) . Ion cyclotron damping will certainly occur if the tur- 
bulent cascade reaches the sm all spatial scales involve d. It 
occurs at the ion inertial scale l|Coles fc Harmonlll989l ). 

Li = 684 km/ v , n E (cm- 3 ) (8) 

and has been clearly observed in the solar wind. Expected 



scales in the ISM range from 300 to 3000 km and it has al- 
most certainly been observed at scales of 300 to 80 km using 
pulse broadening observations l|Bhat et al.l l2004). Damping 
due to ion-neutral collisions is also a resonant process that 
occurs near the ion-neutral collision frequency and would 
result in scales of ^30 AU in typical warm ISM conditions 
l|Minter fc Spangler![l997l ). 



3 OBSERVATIONS AND DATA ANALYSIS 

The PPTA, which commenced observations in Febru- 
ary 2004, uses the Parkes 64-m radio telescope to 
make timing observations of 20 millisecond pulsars. One, 
PSR J1824-2452, lies within the globular cluster M28, the 
others within the disk of our Galaxy. Table [1] lists basic pa- 
rameters for the 20 PPTA pulsars: J2000 name, period (P), 
period derivative (P), orbital period (p,) if the pulsar is in 
a binary system, D M, distance based on the NE2001 elec- 
tron density model l)Cordes fc Lazioll2002h unless the annual 
parallax or another independent distance estimate is avail- 
able, total proper motion (fi) and ecliptic latitude (&e). Be- 
cause PSRs J1022+1001 and J1730-2304 lie very close to 
the ecliptic plane, timing methods cannot provide a precise 
proper motion in ecliptic latitude; for these two pulsars the 
proper motion in ecliptic longitude is given. Each pulsar is 
typically observed at intervals of 2-3 weeks at frequencies 
close to 685 MHz (50 cm), 1400 MHz (20 cm) and 3100 MHz 
(10 cm), where the band designations (based on wavelength) 
are given in parentheses. We have used three backend sys- 
tems: a wideband correlator (WBC), a digital filterbank sys- 
tem (DFBlfl and the Caltech-Parkes-Swinburne Recorder 
2 (CPSR2), a coherent dedispersing system, all of which 
record orthogonal linear polarisations. The WBC provides 
2-bit sampling with a bandwidth of up to 1024 MHz for the 
earlier data. The DFB1 was installed in 2005 June and allows 
8-bit sampling of a 256 MHz bandwidth at 10 cm and 20 cm. 
Observations at 10 and 50 cm are obtained simultaneously 
using a dual-band receiver providing bandwidths of 64 MHz 
at 50 cm and 1024 MHz at 10 cm. Most observations at 20 cm 
are made using the central beam of the Parkes multibeam 
receiver although the "H-OH" receiver has occasionally been 
used for observations in this band. Data are simultaneously 
recorded using the CPSR2 baseband recording system with 
2-bit sampling of two 64-MHz bands, centred on 1341 and 
1405 MHz respectively, and either the WBC or DFB1 with 
256 MHz bandwidth. At 50 cm, data are recorded using one 
band of CPSR2. For all receivers, a linearly polarised broad- 
band calibration signal can be injected into the feed at 45° 
to the two signal probes. 

Observation times per pulsar are typically either 32 min 
or 64 min and data are folded on-line with sub-integration 
times of 1 min for the WBC and DFB1 and 16 s for CPSR2. 
All pulsar observations are preceded by a short (2 min) 
observation of a pulsed calibration signal. For most pul- 
sars the WBC and DFB1 data are split into 512 frequency 
channels with between 256 and 1024 phase bins across the 
pulse period. For CPSR2, the data are coherently dedis- 
persed in each of 128 frequency channels with 1024 phase 

2 A new digital filterbank system (DFB2) with a wider band- 
width and improved resolution is currently under construction. 



4 X. P. You et al 

Table 1. Parameters for the PPTA pulsars. 



PSR P P P b DM Dist. fi b E 

(ms) (lCT 20 ) (d) (cm" 3 pc) (kpc) (masyr" 1 ) (°) 



J0437-4715 


5.757 


5.73 


5.74 


2.6 


0.16 a 


140.89 


-67.87 


J0613-0200 


3.062 


0.96 


1.20 


38.8 


1.71 


7.3 


-25.41 


J0711-6830 


5.491 


1.49 




18.4 


0.86 


21.9 


-82.89 


J1022+1001 


16.453 


4.33 


7.81 


10.3 


0.30 a 


17 


-0.064 


J1024-0719 


5.162 


1.85 




6.5 


0.39 


81 


-16.04 


J1045-4509 


7.474 


1.77 


4.08 


58.1 


1.96 


7.8 


-47.71 


J1600-3053 


3.598 


0.95 


14.35 


52.3 


1.63 


4.1 


-10.07 


J1603-7202 


14.842 


1.56 


6.31 


38.1 


1.17 


8.5 


-49.96 


J1643-1224 


4.622 


1.85 


147.02 


62.4 


2.41 


9 


9.78 


J1713+0747 


4.570 


0.85 


67.83 


16.0 


1.12 a 


6.4 


30.70 


J1730-2304 


8.123 


2.02 




9.6 


0.53 


20.5 


0.19 


J1732-5049 


5.313 


1.38 


5.26 


56.8 


1.41 a 




-27.49 


J1744-1134 


4.075 


0.89 




3.1 


0.36 


20.99 


11.81 


J1824-2452 


3.054 


162.00 




119.9 


3.09 


4.7 


-1.55 


J1857+0943 


5.362 


1.78 


12.33 


13.3 


0.91" 


6.16 


32.32 


J1909-3744 


2.947 


1.40 


1.53 


10.4 


1.14 a 


36.99 


-15.16 


J1939+2134 


1.558 


10.50 




71.0 


3.57 


0.80 


42.30 


J2124-3358 


4.931 


2.05 




4.6 


0.27 


49.0 


-17.82 


J2129-5721 


3.726 


2.07 


6.63 


31.9 


1.36 


8 


-39.90 


J2145-0750 


16.052 


2.98 


6.84 


9.0 


0.50° 


14.1 


5.31 



a Distance obtained from a parallax measurement. 



bins. Off-line proce ssing uses the psrchive software system 
ijHotan et alj liooih . For all recording systems, data from 
frequency channels or sub-integrations which are obviously 
affected by radio-frequency interference are excised, as are 
channels from the outer edges of the band (typically about 5 
per cent of the band at each edge) where the system gain is 
low. Data are then calibrated for variations of instrumental 
gain and phase across the band using observations of the 
pulsed calibration signal and the Stokes parameters formed. 

For all pulsars except PSR J0437-4715, pulse TOAs 
were obtained by cross-correlating a template profile with 
the Stokes I mean pulse profile for each observation. For 
most of the observed pulsars, errors in the calibration pro- 
cedure resulted in TOA errors which were less than the 
uncertainty due to random (receiver) noise. However for 
PSR J0437— 4715 at 20 cm and 50 cm, this was not the case 
and it was advantageous to use the polarimetr ic invariant 
interval instead of Stokes I (|Britton et alJl200Ch . Template 
profiles were formed for each instrument and each observ- 
ing frequency (685, 1341, 1405 MHz for CPSR2, 1369 and 
3100 MHz for the DFB1, 1433 and 3100 MHz for the WBC) 
by weighted summing of all available data to form 'grand 
average' profiles and then blanking the baseline regions. Fig- 
ure [1] shows the grand average Stokes I profiles at the three 
frequencies for all 20 pulsars (except for PSR J2129— 5721 
where we have data for two frequencies only). 

For DM measurements, profile alignment across fre- 
quencies is an important issue. The template profiles for 
a given pulsar were approximately aligned using the cross- 
correlation of each profile with a reference profile. However, 
because of frequency-dependent profile variations, there re- 
mains some uncertainty in the true alignment. In this work, 
we are primarily concerned with variations in DM, not abso- 
lute values, so arbitrary phase offsets between data obtained 
using different systems were included in the timing model. 



The resulting TOAs were analysed using tempo2. Tim- 
ing model parameters were obtained by fitting standard pul- 
sar timing parameters (including astrometric, spin and bi- 
nary parameters) to the 20 cm and /or the 10 cm observa- 
tion^. As we are interested in DM variations, we do not 
fit for any time derivatives of the DM as part of the tim- 
ing model; however, we do allow tempo2 to model the DM 
variations due to the solar wind (this is further discussed in 
H5.2|l . The timing model parameters were subsequently held 
constant in order to obtain timing residuals at all observing 
frequencies. 

To obtain the time variations in DM, ADM(t), we fit- 
ted Equation [5] to segments of timing data. The segment 
lengths were adjusted so that each segment contained at 
least one observation at each frequency (typically 1 or 2 
weeks) . Initially we used the 10 and 50 cm observations to 
obtain ADM(t) because these were obtained simultaneously 
and are well separated in frequency. However, if the pulse 
profile at 10 cm or 50 cm has a low signal-to-noise ratio, we 
also used 20 and 50 cm or 10 and 20 cm to obtain ADM(t). 

The structure function is a useful statistic for studying 
the physical process causing the DM variations. Calcula- 
tion of the structure function is straightforward, even for 
unequally-spaced data (Equation [6]). However, the estima- 
tion of the errors in the structure function due to uncertain- 
ties in the ADM(i) values, and due to the finite duration 
of the ADM(i) series has not been discussed consistently in 
the literature. Because of the irregular data-sampling, r rep- 
resents a "bin" with a width which we have adjusted to give 



3 We use residuals obtained using CPSR2, WBC and DFB1. 
However, for a few pulsars with high DM and short period, 
the WBC profile is significantly smeared and therefore, for 
PSRs J0613-0200, J1600-3053, J1824-2452, J1939+2134, we 
only use the CPSR2 and DFB1 observations. 
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Figure 1. Stokes I profiles for the 20 millisecond pulsars in our sample. For each pulsar we give the 50 cm, 20 cm and 10 cm profiles from 
top to bottom respectively (except for PSR J2129— 5721 where only 50 cm and 20cm profiles are available). 



roughly equal logarithmic sampling. Estimation of Ddm(t) 
for a power-law process from a single "realisation" of the pro- 
cess incurs a signifi cant error which has been discussed by 
iRickett et al.l ([2000 ). For Kolmogorov processes they found 
that the estimation error o cs t (t) oc D(t)[t/(T — t)] 1 ^ s where 
T represents the data-span. We have extended their simula- 
tions to pure power-law processes which do not have a de- 
fined low frequency limit (or "outer scale"). For uniformly 
sampled data from a Kolmogorov process we find that 

<7est(r) = 1.66D(t)(t/T) 1/3 . (9) 

Assuming that the process itself has Gaussian differences, 
the structure function estimator must be ^-distributed, 
and defined by iVdof, the number of degrees of freedom. 
The estimation error can be written in terms of iVdof as 



ff«t(r) = D( T )(2/N doi )° 



N doi = 0.72(T/r) 2 / 3 . As 



data are irregularly sampled we sometimes have fewer pairs, 
N p , or fewer samples, N 3 , contributing than A^of ■ We there- 
fore approximate the actual number of degrees of freedom 
as the minimum of (N P ,N 3 , Ndoi)- 

In addition to the estimation error discussed above, 
which is the amount by which a single realisation of a ran- 
dom process can be expected to depart from the theoretical 
mean, we must consider the fact that the measurements in- 
clude a white noise component independent of DM(t). We 
assume that the errors on the ADM(t) values are indepen- 
dent, Gaussian and have known, but different, standard de- 



viations. These contribute a different bias and error to each 
D est (T) which are computed by expanding each D ea t(r) as 
shown in Appendix A. 

Our work contrasts with earlier estimates of the struc- 
ture function. In earlier work, the uncertainty on the struc- 
ture function was either estimated as the standard devia- 
tion divided by th e square root of the number of points 
l|Cordes et al.lll990l) or as the smaller of the number of points 
and {T/t) 1 Ramachandran et al.l [2006 s ) . The choice of the 
bias term which is subtracted from Ddm (t) has also been in- 
consistent in the literature. Usually the bias has been taken 
as the average of the ADM(i) errors. This method is only 
accurate when the number of points is very large and the 
error on ADM(i) is significantly smaller than its value. Ear- 
lier work also did not allow for the error on the measured 
ADM(i) values which, for some data-sets, is very important. 

It is useful to obtain the diffractive time-scale, Td, and 
bandwidth, Vd, for each pulsar to compare with the -Ddm(t). 
For many pulsars we were able to obtain these values from 
the literature. Table [2] gives these Td and Vd values. In col- 
umn order, the Table contains the pulsar name, observing 
frequency, Td, Vd and a bibliographic reference. However, no 
measurements existed for nine of our pulsars. As our data 
have relatively poor frequency and time resolution (for this 
purpose) , we obtained r d an d v d using a structure function 
analysis ( (Rickett et alj|2000| ) rather than the more standard 
method of using the auto-correlation function of the dy- 
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Table 2. Scintillation parameters for the PPTA millisecond pul- 
sars 



PSR Name 


Freq. 
^IViriz 1 


Td 
(min) 


Vd 

(IVlriz ) 


Iter. 


J0437- 


-4715 


327 


1.9 - 5.1 


0.18 - 3.0 


i 






436 


4.6 - 11 


3.2 - 4.4 


2,3 






660 


7.8 


17 


2 


J0711- 


-6830 


436 


13 


0.37 


2 






660 


16 


1.2 


2 


J1600- 


-3053 


1373 


4.7 


< 0.5 


4 


J1603- 


-7202 


660 


9.2 


0.36 


2 


J1713- 


-0747 


430 


14 


0.6 


5 






436 


28 


1.5 


2 


J1730- 


-2304 


327 


7.4 - 7.5 


0.10 - 0.12 


1 






436 


6.3 - 12 


0.15 - 0.18 


2,3 






660 


9.7 


1.4 


2 






1520 


23-27 


30 - 38 


2,3 


J1744- 


-1134 


436 


21 


1.3 


2 






660 


20 


2.3 


2 


J19394 


-2134 


320 


1.1 


0.0014 


6 






430 


1.7 


0.0042 


6 






1400 


7.4 


0.92 


6 


J2124- 


-3358 


436 


44 


6.9 


2 


J2129- 


-5718 


436 


11 


0.29 


2 






660 


17 


1.3 


2 






1520 


24 


58 


2 


J2145- 


-0750 


327 


6.4 


0.33 


1 






436 


21-25 


0.61 - 2.5 


2,3 



Refere nce: (1) iGothoskar fc Guptal d2000f); (21 I Johnston et al.l 
lll998ll: (3llNicast.ro fc Johnston! lll995ll: (4) IOrd et al.l j2006t l: C5) 
iBogdanov et alj J2002t ); (6) ICordes et alj Jl990l) . 

namic spectrum. The implementation of this technique is 
discussed in Appendix B. From each Td measurement wc 
have used the definition that D^(jd) — 1 to obtain an esti- 
mate of 

D DM {T d ) = {Ku/2tv) 2 . (10) 



4 RESULTS 

The variation of DM with time, ADM(t), for each pulsar is 
shown in Figure [2] We list, in Table [3] the bands used for 
each pulsar and the interval over which the ADM(t) values 
were measured and the slope of the best-fitting straight line 
dDM/di. The panels in Figure [2] are chosen so that all pul- 
sars have the same time axis, but different scalings are used 
for the ordinate. We see large-scale DM variations for six of 
our pulsars (PSRs J0437-4715, J1045-4509, J1643-1224, 
J1824-2452, J1909-3744, J1939+2134) with a maximum 
range in ADM of ~ 0.014 cm" 3 pc for PSR J1045-4509. 

4.1 Diffractive scintillation parameters 

Using the method described in §3, we obtained diffractive 
scintillation timescales for 17 of our pulsars, obtaining Td 
and Vd values for each observation with a high S/N. In col- 
umn order, Table U contains the pulsar name, observing fre- 
quency and the range of our measured Td and Vd values, 
respectively. For pulsars where previous measurements exist 



Table 3. Summary of the DM variations for the PPTA pulsars 



PSR Name 


Band 


Interval 


dDM/dt 




Data Span 




(cm) 


(d) 


(cm _J pc yr — 




(yr) 


J0437— 4715 


10, 


50 


15 


1.0(2) 


X 


10 




3.0 


JUolo— U2UU 


20. 


50 


15 


-9(2) 


X 


10 




2.9 


JlU 11 — OOoU 


20. 


50 


15 


-2.5(9) 


X 


10 


g 


2.8 


J lUzz+lUUl 


10. 


50 


7 


-5(60) 


X 


10 


-6 


2. t 


T i no a r\'7~\ n 


20. 


50 


15 


3.2(9) 


X 


10 


-4 


2.9 


J1U45— 45U9 


20. 


50 


15 


-5.56(9) 


X 


10 


-3 


2.9 


J 1DUU — OUOo 


10. 


20 


10 


-9(2) 


X 


10 


-4 


9 7 


J1603-7202 


20, 


50 


15 


1.28(5) 


X 


10 


-3 


2.6 


J1643-1224 


20. 


50 


15 


-1.18(6) 


X 


10 


-3 


2.6 


J1713+0747 


20. 


50 


15 


5(23) 


X 


10 


-6 


2.7 


J1730-2304 


20. 


50 


7 


3.9(8) 


X 


10 


-4 


2.5 


J1732-5049 


20. 


50 


7 


-7(1) 


X 


10 


-4 


2.4 


J1744-1134 


20. 


50 


15 


5(2) 


X 


10 


-5 


2.7 


J1824-2452 


20. 


50 


12 


6.0(1) 


X 


10 


-3 


1.1 


J1857+0943 


20. 


50 


15 


1.5(1) 


X 


10 


-3 


2.2 


J1909-3744 


20. 


50 


15 


-3.28(6) 


X 


10 


-4 


2.7 


J1939+2134 


20. 


50 


15 


2.57(2) 


X 


10 


-4 


2.3 


J2124-3358 


20. 


50 


15 


2.5(8) 


X 


10 


-4 


2.7 


J2129-5721 


20. 


50 


7 


-2(3) 


X 


10 


-4 


0.9 


J2145-0750 


20. 


50 


15 


4.0(4) 


X 


10 


-4 


2.1 



our results are consistent with values in the literature. The 
scatter in Td observations is much greater than the error 
bars on individual Td measurements. We have confirmed, by 
simulation, that the reason for this is that Td is estimated 
from observations which are much shorter than the refrac- 
tive scale. 

For three pulsars it was difficult for us to obtain 
diffractive time-scale measurements. For PSR J0437— 4715, 
this is not a problem as there are many measurements 
available in the literature. The diffractive time-scale for 
PSR J1824-2452 is too short at 20 cm (less than 1 min) for 
us to measure. At 10 cm, this pulsar is very weak (SNR ~20 
for a 1 hr observation), but we were able to obtain a few us- 
able observations. The Td for PSR J2124— 3 358 is relatively 
long. From the Td = 44 min at 436 MHz lijohnston et al.l 
Il99ct l. we can estimate that Td at 685 MHz is ~76min, but 
our current observation time for this pulsar is only 32 min. 

4.2 Structure functions 

We have calculated structure functions from ADM(t) for 
each of our pulsars. Representative examples are shown in 
Figure [3] and Figure [4] In these figures, we have included Td 
measurements obtained from the literature (cross-symbols) 
or from our data (open triangle symbols). For some pul- 
sars, we have been able to derive an estimate of Ddm at 
large time lags from dDM/dt measurements in the literature 
( Ho bbs et ai1l2004h that were obtained using a single data- 
selQ; such points are indicated using a full square-symbol at 
the rightmost edge of the plot. To put the data in context, we 
have drawn two theoretical lines fitted through the Td points, 

4 A given dDM/dt measurement can be converted to a single 
point on a structure function as (dDM/dt ■ T) 2 , where T is the 
data span. 



DM Variations and Pulsar Timing 7 



CO 

o 

CO 

I 

03 
CO 

O 
03 
CO 

I 

03 
CM 

O - 

03 
CM 
I 



CO 
CM 



CO 
CM 
I 



O h 

in 



b I 



Q 



CM 
CO 



o - 

CM 
CO 
I 

CM 
O - 

CM 
I 



CO 

o 



CO 

I 



o 

I 

CO 

o 

CO 

I 



-1 1 1 1 1 1 1 1 1- 



H — ' — ' — ' — ' — I — >- 



J0437-4715 



J0613-0200 

H — 1 — 1 — 1 — ' — I — >- 



♦♦1 

H 

I I I 



— f I T l J1730-2304 " 

t 1 — "-r 1 — ' — 1 — I — 1 — 1 — ' — 1 — I — 

J1732-5049 J " 

H — i — i — i — i — 1—i — i — j — i — I — 



J071 1-6830 

1 — i — i — i — i — I — "r 



i i I 



J1022+1001 



— | 1 1 1 1 1 1 1 1 1 1- 

L9 



J1024-0711 
I i i i i I 



J1045-4509 



i i I 




J1603-7202 
-| 1 1 1 1 1 1 1 1 1 1- 



1 v v * A 

J1643-1224 * 
H — 1 — 1 — r 1 — 1 — h — 1 — 1 — 1 — 1 — f 



♦ 



— i — i — p — i — r — 1 — 1 — 1 — 1 — I — 
"17 " 



J1713+0747 
i i i i I u 



-i 1 1 1 ill 1 1 r 



J1744-1134 
I 1 1 1 1 I 



i i i I 



- CM 
I 



J1824-2452 
-| 1 1 1 1 1 1 1 1 1 1- 



.» ♦ . . ' 

** * J1857+0943 * 
-| 1 1 1 1 1 1 1 1 1 1- 



i* . 

J1909-3744 T +*.% H 
-| 1 1 1 1 1 1 1 1 1 \ 



J1939+2134 
-| 1 1 1 1 1 1 h- 



+ 



J2 124-3358 
-| 1 1 1 1 1 1 1 1 1 f 



V 



J2129-5721 
H 1 1 1 1 1 1 \ 1 — ; — H — jf 



t - - 

♦ -I o 



J2145-0750 



o 

CM 

iH 

I 

T-H 

o 

i-i 
I 

CM 
O 



CO 



- o 



- CO 

- £ 

CO 00 

I 

Ho 6 

o 

-I CO 7 
I o 

CO 



- o 

CO 

I 

CO 
CO 



- o 

CO 

-\ CO 

I 

CM 
O 



- CM 
I 

03 
CO 

- O 
03 

- CO 
I 



-500 



500 -500 



500 



MJD-53500 



MJD-53500 



Figure 2. DM variations of 20 millisecond pulsars. Note a ADM of 10 4 cm 3 pc corresponds to time delays at 10 cm of 43 ns, at 20 cm 
of 212 ns and at 50 cm of 884 ns. 
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Figure 3. Structure functions, Ddm(t) for four pulsars. The derived estimates are marked by triangles from our data and crosses 
from the literature. The estimates obtained directly from the time series ADM(t) as discussed in Appendix A are marked as filled circles 
with error bars. Open circles indicate a negative estimate. A Kolmogorov model fit to r d is shown using a solid line and a quadratic 
model is shown dash-dotted. Confidence limits on the Kolmogorov model are solid lines bracketting the model. Equivalent delays at 
1400 MHz are shown for 100 ns and 1 fj,s. For PSR J1022+1001, a point derived from dDM/dt is shown as a solid box with error bars at 
the longest time lag. 



one (full line) with the Kolmogorov exponent (a = 5/3), 
the other (dashed-dotted line) with a = 2 corresponding 
to the steepest po ssible structur e function resulting from 
plasma turbulence (| Rickettl 1 1 9 9 Oh (hereafter, this spectrum 
is known as "quadratic"). Estimation error bounds on the 
theoretical Kolmogorov model (at the 68% confidence level) 
for each data point are plotted using solid lines that bracket 
the theoretical model. 

The remaining symbols used on the figures are as fol- 
lows. The structure function values measured from ADM(i) 
are plotted using solid circle symbols. The errors on these 
points are estimated from the uncertainty on ADM(i). For 
cases where the error is larger than the value we use down- 
ward pointing arrow symbols for the lower bound on the 
error bar. As we subtracted the bias due to the uncertain- 
ties on the ADM(i) measurements, it is possible for large 
uncertainties on ADM(i) that the structure function values 
are negative. We indicate such points using open circles and 
a downward arrow plotted at Ddm(t) plus twice its error. 
The structure function plots all have the same scaling. 

For comparison, we also indicate the value of Ddm 
that would be expected for white timing residuals with a 
given rms (<r rms ) at 1400 MHz. The relationship between 



the structure function of the timing residuals, Dtoa, and 
the structure function of DM variations, Ddm, is 



-Ddm = Dtoa(Ku 



(11) 



If the timing residuals are white, then Dtoa = 2a^ m3 . We 
indicate, using solid horizontal lines, white noise with an rms 
of 1 [is and 100 ns. 

For PSRs J1045-4509, J1824-2452 and J1909-3744 
we have added an indication of the inner time-scale (see 
§5.3). For PSRs J1939+2134 and J1824-2452, we also over- 
lay a dotted line which is the structure function from earlier 
work (see §5.3). 

Our Td valu es for PSR J0437— 471 5 can be compared 
with the wo rk of |Smirnova et al. | (|2006l )_ who scaled the ob- 
servations of lJohnston et aL I (|l998l ) and lGothoskar fc Gupta! 
(2000) by a large factor, assuming that the scintillation 
index was much smaller than unity, although the scintil- 
lation index w a s not reported by the original observers. 
ISmirnova et al.l l|2006h deduced a phase structure function 
which is two orders of magnitude lower than ours in the 
vicinity of Td. 
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Figure 4. As for Figure|3] but for another four pulsars. For PSRs J1643— 1224 and J1939+2134, a point derived from dDM/dt is shown 
as a solid box with error bars at the longest time lag. 
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Figure 5. DM variations of PSR J1022+1001. The dashed curve is the TEMP02 modelled DMq values. Panel a shows the DM variations 
with no correction for the solar wind. Panel b shows the DM variations after correction by the model used in TEMP02. 



4.3 Summary of results 

For all the pulsars we find that the measured structure 
functions lie above the lower error bound of the Kol- 
mogorov model. Two, PSRs J0437-4715 and J1939+2134, 
fit the Kolmogorov model well. Two, PSRs J1045-4509 



and J1909— 3744, are clearly inconsistent with a pure Kol- 
mogorov power law, requiring a large inner scale. One, 
PSR J1824— 2452, has few good Td measurements, but may 
well have Ti > t&. The remaining 15 pulsars are dominated 
by white noise at small lags, and for five of these we can- 
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Table 4. Scintillation parameters from our observations 



PSR Name Freq. 

(MHz) (min) (MHz) 



J0613- 


-0200 


1369 


10 




54 


0.98 




3.1 


J0711- 


-6830 


685 


18 




42 


1.0 




5.4 






1369 


36 




127 


30 




77 


J1022-1 


-1001 


685 


53 




18 1 


6. i 




10 


J1024- 


-0719 


685 


23 




89 


4.8 




17 


J1045- 


-4509 


3100 


2.2 




12 


0.64 




15 


J1600- 


-3053 


3100 


1.0 




23 


1.90 




6.9 


J1603- 


-7202 


685 


8.6 


— 


27 


1.5 


— 


7.8 






1369 


7.7 




40 


1.8 




18 


J1643- 


-1224 


3100 


2.6 




12 


0.89 




2.0 


J1713+0747 


685 


20 




47 


1.8 




11 


J1730- 


-2304 


685 


9.3 




28 


1.2 




5.4 






1369 


17 




54 


3.9 




32 


J1732- 


-5049 


1369 


18 




39 


1.8 




6.2 


J1744- 


-1134 


685 


29 




76 


1.1 




40 


J1824- 


-2452 


3100 


1.1 




9.5 


0.6 




1.1 


J18574 


-0943 


685 


13 




22 


2.5 




7.8 






1369 


16 




68 


2.7 




25 


J1909- 


-3744 


685 


16 




69 


2.8 




20 


J19394 


-2134 


1369 


4.1 




10 


1.8 




5.4 


J2129- 


-5718 


685 


12 




39 


1.7 




4.5 






1369 


35 




79 


25 




234 


J2145- 


-0750 


685 


20 




111 


2.9 




44 



not constrain the slope of the underlying power-law spec- 
trum. Two pulsars (PSRs J1744-1134, J1857+0943) could 
not be classified on the basis of our measurements, but 
appear to be Kolmogorov on the basis of previously pub- 
lished dDM/dt values. For five pulsars (PSRs J0613-0200, 
J1600-3053, J1643-1224, J1713+0747 and J1732-5049) 
the structure functions fall below the quadratic model 
at large time lags, strongly suggesting that the under- 
lying spectrum is Kolmogorov. The final three pulsars 
(PSRs J1603-7202, J1730-2304 and J2124-3358) appear 
to follow the quadratic model at large lags. However it 
should be realised that the structure functions at large lags 
are relatively poorly estimated and this separation of the 
pulsars into different categories is not perfectly clear. 

The results described above lead us to propose that the 
structure functions for all our pulsars contain an ISM com- 
ponent that is either a pure Kolmogorov power-law or a 
Kolmogorov power-law with a large inner scale. 



5 DISCUSSION OF DM MEASUREMENTS 

5.1 Comparison with earlier work 

Much of the earlier work has concent rated on measuring 
(and modelling) d DM/dt values (e.g. iBacker et all 1 19931 ; 
iHobbs et al. I l2004l) . For comparison with earlier work we 
have listed in Table [3] the slope of the best-fitting straight 
line across the entire data-set for each of our pulsars, 
dDM/dt. Our values generally do not agree with the pre- 
viously published values. However, for our data-sets, a sin- 
gle dDM/dt value models the observed ADM(i) values well 
only for a few pulsars (PSRs J1045-4509, J1824-2452 
and J1909-3744) and the dDM/dt values for other pulsars 



are misleading. For instance, for PSR J0437— 4715 our re- 
sults indicate that the DM evolution for this pulsar can 
roughly be described using three dDM/dt values: prior to 
MJD 53400 dDM/dt = (-2.98 ± 0.07) x 10 _4 cm~ 3 pcyr _1 , 
between MJD 53400 and 53700, dDM/dt = (6.2 ± 0.2) x 
10~ 4 cm~ 3 pcyr _1 and subsequently dDM/dt = (3 ± 1) X 
10 _5 cm _3 pcyr _1 . Clearly, the ADM(t) values are better de- 
scribed using the structure function. 

5.2 The solar wind 

The solar wind leads to a significant change in DM for pul- 
sars close to the ecliptic plane during close approaches of the 
line of sight to the pulsar with the Sun. According to the 
tempo2 model described in §2, an unmodelled solar wind 
contributes ~ 100 ns at 1400 MHz for sources within 60° of 
the Sun and ~1 fis within 7°. It is therefore clear that correc- 
tions are necessary for 18 out of our 20 PPTA pulsars. The 
correction can potentially be made by modelling the solar 
wind or by directly measuring the DM to sufficient accuracy 
using multiple frequency observations. 

The solar wind varies with time and position. An 
overview of th e relevant solar physics can be found in 
ISchwennl (|2006l ). Most of the variations in the solar wind 
are ascribed to a slowly changing spatial pattern that ro- 
tates with a 27-day period. In addition global transients, 
called coronal mass ejections (CME), occur every few days. 
The chances of a given observation, which typically has 
a duration of 30-60 min, encountering a CME are only a 
few percent, so these are not the most important effects. 
The quasi-static spatial pattern is roughly bimodal, the 
"slow solar wind" w ith high electron densi ty is concentrated 
within about ±20° (|McComas et al.ll2000l ) and the "fast so- 
lar wind" with lower electron densities at higher latitudes. 
The density difference is a factor of four at 1 AU and in- 
creases near the Sun. The correction necessary for a given 
observation can then vary by a factor of four depending on 
how much of the line of sight is in the slow versus the fast 
wind. 

Corrections for the solar wind have been attempted 
in both tempo and tempo2. Both programs implement 
constant-density spherically symmetric models. Tempo uses 
a high density model where the electron density at 1 AU, 
n e (lAU) = 10cm -3 , whereas t empo2 has a lower den- 
sity model of n e (lAU ) = 4cm~ 3 . ISplaver et all (|2005l ) and 
lLommen et al.l (|2006l ) used an identical spherically symmet- 
ric model, but fitted for the electron density. They obtained 
ne(lAU) = 5±4cm" 3 and n e (lAU) = 6.9±2.1cm~ 3 respec- 
tively. However, it is not possible for a spherically symmetric 
model to correct the average timing residual due to the large 
difference in density between the fast and slow winds. This 
is clearly demonstrated by our PSR J1022+1001 observa- 
tions. The left panel in Figure [5] shows the DM variations of 
this pulsar without correcting for the solar wind and gives 
the correction from the tempo2 model as a dashed line. 
The right panel shows the DM variations after correction 
using the tempo2 model. During the year 2004, the tempo2 
model did accurately correct the effect. The original tempo 
model which uses a larger electron density overcorrects these 
observations. The opposite occurs during 2005, when the 
tempo2 model under-corrects the observations whereas the 
original tempo models the solar wind well. 
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It is possible to use coronal measurements to improve 
our correction by estimating which parts of the line of sight 
are in the fast and which in the slow wind. This can be 
demonstrated with our PSR J1022+1001 data, but it is not 
yet clear whether using an updated model to correct the 
observations improves on simply measuring the excess DM 
using mult i- frequency observations. This work will be re- 
ported in a future publication. 



5.3 Spectrum of the ISM 

All but six pulsars are consistent with a Kolmogorov fluctua- 
tion spectrum with an inner time-scale smaller than Td- The 
clearest examples are PSRs J0437-4715 and J1939+2134. 
The structure function for PSR J0437-4715 lies slightly 
above the upper bound of the Kolmogorov model fit to the 
Td data. However, if the Td data were divided by a factor of 
1.35 they would be consistent. A line with this shift is shown 
in Figure [3^ through the Td data. Given the large scatter in 
Td we consider that there is no evidence for an inner scale, 
nor do we see a need to rescale D^ijd) as did Smirnova et 
al. (2006). The structure function at the largest time-scales 
is currently consistent with a Kolmogorov process, but there 
is an indication that the structure function may be flatten- 
ing at these scales, as if an outer scale around 60 AU were 
present. We do not expect such a small outer scale, but it 
is not impossible if the turbulence has inhomogeneities of 
this order. Such structures could be caused, for example, by 
shear instabilities or due to a large-scale damping mecha- 
nism such as ion-neutral damping. The presence or absence 
of this flattening will become clearer in a few years, as we 
accumulate more observations of this pulsar. 

There have been several analyses o f the struc- 
ture function for PSR J19 39+2134 (jCordes et all 
1990l; iRamachandran et alJ [2006). The recent result of 
Ramachandran et al.l (120061 ) is overlaid on our structure 
function in Figure [4}i. They fitted for the power-law 
exponent and obtained a = 1.66 ± 0.04. They also compare 
their Ddm with a single Td measurement of 180 s. This 
comparison suggests an inner scale of 1.3 x 10 9 m. Our Ddm 
observations are slightly above the Kolmogorov model fit 
to the Td observations. However, as with PSR J0437— 4715 
there is a large scatter in Td- A shift of 1.43, shown as a 
solid line through the data in Figure [3JI, would make the 
Ddm consistent. Thus we believe that the case for an inner 
time-scale greater than Td is weak for this pulsar. 

The structure functions for PSRs J 1045— 4509 and 
J1909— 3744 both lie well above the upper bound of the Kol- 
mogorov spectrum and require an inner scale which is much 
larger than Ut^. We can identify a break in the structure 
function for PSR J1045— 4509 which suggests an inner time- 
scale of about 12 d. The observations of PSR J1909-3744 
can set a lower bound on the inner time-scale of 500 d. To 
our knowledge these are the first observations of such large 
inner scales. In order to estimate the corresponding spatial 
scales we use V » 3 .85 x 10 4 (i/dD) 1 ^ 2 / (vTd) for a thin screen 
i|Gupta et al.lll994h . In both cases the inner scales of 0.7 and 
20 AU are comparable to or larger than the refractive scales 



of 0.38 and 0.19 AU respectively Q. These scales are so large 
that they can only be compared with ion-neutral damping 
scales. Ion-cyclotron damping at such scales would require 
absurdly low densities (see Equation [8} . 

Measurements of the structure function 
for PSR J1824- 2 452 have been published by 
ICognard &: Lestradei l|l997h and shown to be consis- 
tent with a Kolmogorov process. In Figure |4p we overlay 
the earlier structure function (dotted line) with our results. 
The two analyses are consistent. Our Td estimates suggest 
that n > Td, but because there are few measurements of Td 
and the inferred n is not as a large as for PSRs J1045— 4509 
and J1909— 3744, the evidence for a large inner scale is 
weak. 

In analysing the structure functions we have assumed 
that the line-of-sight velocity is constant. However the true 
velocity is a vector sum, weighted by the distance of the 
scattering plasma, of the pulsar proper motion, the orbital 
velocity for binary pulsars, the velocity of the plasma with 
respect to th e local standard of r est and the orbital velocity 
of the earth l|Rickett et al.ll2000h . For many of our pulsars 
these effects are important. Thirteen of our pulsars are in 
binary systems, and in five of those the orbital velocity is 
comparable with the proper motion. In all five of these plus 
another four non-binary pulsars the Earth's orbital velocity 
is also comparable. For these pulsars the magnitude and 
direction of the velocity can change significantly, both on a 
time-scale of days, and annually. 

Diffractive observations are made on a time-scale which 
is short compared to the orbital periods, so such observations 
are affected by the instantaneous vector sum of velocities. 
DM variations are normally averaged over times longer than 
the typical binary periods in our sample, but shorter than a 
year. Thus these measurements are not affected by the bi- 
nary orbital velocity. On average the diffractive observations 
see a higher velocity than the DM observations which will 
lead to the temporal structure function being flatter than 
the spatial structure function. If the turbulence in the inter- 
stellar plasma is anisotropic, and there is increasing evidence 
that such anisotropy is common, then the apparent diffrac- 
tive time-scale will depend strongly on the direction of the 
velocity. 

The net effect on our estimation of the structure 
function of DM is not large, because four of our five 
best constrained observations have relatively small veloc- 
ity modulation. However, observations of the solitary pul- 
sar PSR J1939+2134 are strongly modulated by the Earth's 
orbital velocity. In this case both diffractive and DM ob- 
servations see the same time-varying velocity so the slope 
of the structure function is not altered, but both observa- 
tions will be "noisier" than expected. In fact, the Ddm(t) 
for this pulsar is noisier than expected, suggesting that we 
are seeing the effect of annual variations in velocity. This 
effect is largest in PSR J2145— 0750, for which the proper 
motion and orbital velocity are both similar to the Earth's 
orbital velocity. This pulsar shows a wider spread in Td than 
most. The structure function for this pulsar is white-noise 
dominated, making it difficult to estimate its slope of the 



5 The refract ive time-scales are found using t r 2vlvd T d 
llRickedll990l) 
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structure function. The expected velocity modulation makes 
it even more difficult , so even though the structure function 
appears to be Kolmogorov on the basis of a single dDM/dt 
value, more observations are required to confirm this. 

In PSRs J1045-4509 and J1909-3744 we observe inten- 
sity variations at the diffractive time-scale. However if the 
inner scales were greater than the refractive scales as they 
appear to be, then the structure function is quadratic and n o 
intensity scintillation should be observed* ^ Wandzurall 1980h . 
Thus there must be an underlying Kolmogorov process 
which is roughly equal in amplitude to the steep-spectrum 
process at r<j. 

This situation has been proposed theoretically (Zweibel, 
private communication 2006). It can arise when the primary 
energy input to the turbulence is at scales larger than the 
ion-neutral damping scale. Energy will cascade down to the 
ion-neutral scale where most of it will be absorbed. However 
some energy may 'tunnel' through the damping region to 
support a second Kolmogorov cascade at a lower level (in 
the vicinity of the ion-neutral collision frequency, plasma 
waves are evanescent). 

In the case of PSR J1045— 4509 the energy difference 
is about a factor of 30. For PSR J1909-3744 we can only 
say that the energy difference is at least a factor of 30. This 
is a very intriguing possibility which requires both observa- 
tional and theoretical followup. The existence of steep spec- 
tra might be confirmed by VLBI observations which should 
show an rms position wander of X/(2xVrd) on a time-scale of 
the inner time-scale (where A is the wavelength). The base- 
lines needed for a 1 rad rms phase difference are Vr<j, where 
r d scales as is 12 . This is about 8000 km for PSR J1045-4509 
and 6000 km for PSR J1824-2452 at 1400 MHz. It would 
be much more difficult to measure PSR J1909— 3744 as the 
baseline, even at 327 MHz, would be 50000 km and the time- 
scale for position wander would be greater than 500 d. 

5.4 White noise 

Twelve of our pulsars show a well-defined flattening of the 
structure functions at small lags which indicates the pres- 
ence of a white-noise process that is substantially greater 
than the measurement error. Although this has not been 
discussed in the context of DM measurements before, it is a 
well-known anomaly in TOA measurements. Observers have 
often rescaled their measurement error estimates to match 
this white noise, but it is not clear that the additional white 
noise is due to measurement error. It could be due to a 
process intrinsic to the pulsar, unexplained calibration is- 
sues, or to diffractive TOA noise. The latter process will 
have the same time scale as diffractive intensity scintilla- 
tion and is highly correlated with the intensity scintilla- 
tion. Its rms is of the order of to = 1/(2-7^). Although 
this phenomen a has not been well-s tudied, it has been ob- 
served dfiectIy_^e^^^^t^L|[l99l) and discussed theoret- 
ically (|Romani et al.lll986l ). We have compared the observed 
white noise rms with to for each pulsar and find that this 
mechanism will need to be considered for four of the PPTA 



6 A pure quadratic structure function corresponds to a linear 
phase gradient, which simply shifts the apparent position of the 
source and does not change the intensity. 



pulsars: PSRs J1045-4509, J1600-3053, J1643-1224 and 
J1939+2134. Since this mechanism is correlated with inten- 
sity it may be possible to use intensity measurements to 
correct it. It is likely that most of the white noise is related 
to system calibration errors as it is known, at least for some 
pulsars, to depend on the observing frequency and backend 
instruments. 



6 CORRECTION OF RESIDUALS FOR DM 
VARIATIONS 

During the design of the PPTA project it was realised that 
DM variations would be an important source of timing noise 
and, unless corrected, would obscure the signature of many 
interesting phenomena such as gravitational waves. Initial 
expectations were that observations using the dual-band re- 
ceiver would be used to determine ADM(i) which, in turn, 
would be used to correct the 20 cm timing residuals. Of 
course, the measurements of DM include a white-noise com- 
ponent discussed earlier and hence, the correction for the 
"red" DM variations adds white noise. Smoothing the DM 
data before making the correction will reduce the white noise 
more than the DM noise. However, choosing the optimal 
smoothing is non-trivial. 

The problem is that the timing model includes numer- 
ous terms such as parallax, position, proper motion, period 
and period derivative that absorb some of the residuals due 
to DM variations. Fitting the timing model to the residu- 
als can substantially reduce the effect of DM variations, but 
causes significant errors in the fitted parameters. We cannot 
use the post-fit rms timing residual as a goodness measure, 
since it would not change at all after correction if the DM 
effect has been totally absorbed in the fitted parameters. 
Accordingly we have calculated the optimal smoothing fac- 
tor using a simple analytical model which requires equally 
spaced observations. We have also simulated the observa- 
tions at the actual sampling intervals with known model 
parameters and determined the optimal smoothing for the 
simulated data without having to fit a timing model. This 
work is outlined in Appendix C. 

An example of the correction process is shown in 
Figure [6] We have shown the post-fit residuals for 
PSR J1939+2134 before correction (with an rms timing 
residual of 0.291 fis), and after correction (rms of 0.193 fjs). 
We can assume that the timing model after correction is 
much more accurate than before correction. Therefore we 
have plotted the residuals of the uncorrected data using 
the more accurate corrected model resulting in an rms of 
0.941 [is. The total proper motion, pulse-frequency, and fre- 
quency derivative were changed during the correction pro- 
cess by 7a, 28a and 28cr respectively. 

Figure|6]shows clearly the necessity of correction for DM 
variations, and it also shows how fitting a timing model can 
spuriously remove a "red" process. For instance, if the TOA 
variations included the signature of a gravitational wave 
which resembled the DM variations (both are expected to 
have a steep, "red" signature) then fitting the timing model 
would also have removed most of the gravitational wave sig- 
nature. Fortunately as the duration of the timing data in- 
creases it becomes harder for the timing model to emulate 
either the DM variations or the signature of gravitational 
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Table 5. Improvement of timing residuals after correction for the 
DM variations. 
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Figure 6. Timing residuals for PSR J1939+2134. The upper 
panel (a) shows the timing residuals before DM correction, panel 
(b) contains the timing residuals after correcting for the DM vari- 
ations using a 71 day smoothing. Panel (c) shows the residuals 
obtained using the corrected parameters with the uncorrected 
TOAs. 



waves. This is because the terms related to the motion of 
the Earth have annual or semiannual periods so they are 
not as effective at removing longer period variations. 

The process described above has also been applied to 
five other pulsars for which DM fluctuations are important. 
For each of these pulsars, the theoretical smoothing time for 
uniformly sampled data and the actual optimal smoothing 
time from the simulation are given in the first two columns of 
TableO The rms timing residuals corresponding to the three 
panels of Figure |SJ "original" , "corrected" , and "true un- 
corrected" (where the corrected pulsar parameters are used 
to model the uncorrected TOAs), are tabulated in the last 
three columns of the table. There is a correlation between 
the improvement of the residuals and the slope of ADM(t). 
The pulsars with the least slope, PSRs J0437— 4715 and 
J1939+2134, showed the most improvement in rms after 
correction for the DM variations. This is because a linear 
slope can be corrected exactly, but spuriously, in the origi- 
nal fitting. The ADM(t) values for all except the last pul- 
sar in the table are dominated by the plasma contribution. 
These all show significantly higher "true uncorrected" resid- 
uals demonstrating that the correction process was impor- 
tant even though it may not have significantly lowered the 
rms timing residuals. The last pulsar, PSR J1643— 1224, is 
dominated by white noise and does not show much improve- 
ment in residual, nor is the true uncorrected residual much 
larger. As the PPTA continues to collect more data, the DM 
corrections will become increasingly important and will have 
to be applied to more of the observed pulsars. 



7 CONCLUSIONS 

We have shown that correction for the plasma delay is es- 
sential for the purposes of the PPTA project and have de- 
veloped an optimal way of applying this correction. 

We also show that the spherically symmetric solar wind 
models included in the pulsar timing packages tempo and 
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tempo2 are of marginal value. More sophisticated models 
may be useful, especially in situations where it is difficult or 
impossible to measure the DM variations directly. 

We have analysed the observed DM variations and 
found that most are consistent with a simple Kolmogorov 
model of interstellar turbulence with dissipation at a rela- 
tively small scale such as would be caused by ion cyclotron 
damping. However at least two of the 15 pulsars for which we 
can estimate the spectral exponent require a steeper spec- 
trum and suggest strongly that ion-neutral collisions are im- 
portant in damping the turbulence spectrum at AU scales. 
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APPENDIX A: CALCULATION OF THE 
STRUCTURE FUNCTION 

Let D cst (r) be the estimated structure function for the DM 
variations at time lag r. As described in §2, we assume that 
the errors on the ADM(t) values are independent, Gaussian 
and have known, but different, standard deviations. The 
bias and error that these errors contribute to D C st(i~) are 
obtained by expanding each D C st(i~) as 



D cst (T) 



= w E[ ADM «- ADM ^] 2 

P l ij 

+ ^>(*) 2 + e(jf] 

ij 

+ 2$}e(t) - e(j)][ADM(i) - ADM(j)] 

ij 

ij ) 



(Al) 



where N p is the number of pairs, ADM(i) is the ADM within 
the time lag "bin" r , and e(i) is the corresponding error. 
The first term in this expansion is the desired estimator and 
the other terms are uncorrelated errors. The second term is 
the only error term that does not have zero mean and so 
contributes a bias which must be calculated and subtracted 
from Ddm(t). The variances of each term are easily calcu- 
lated and summed to give the total variance in Ddm(t). So 
finally, the calculated Ddm(t) is 



D 



dm(t ) 



= ^|e[ adm «- adm w] 2 

P I ij 
- X> l£ (i) 2 j 



(A2) 



where e(i) is the rms of the e(i), Ni is the number of times 
that ADM(i) is used to calculate the -Ddm(t). The variance 
of the Ddm(t) is 

+ 4EE e « 2 [ ADM «- ADM W] 2 



+ 4E^) 2 c(.?) 2 } 

ij ) 



(A3) 



We have developed a tempo2 plug-in that is publicly avail- 
able (download and usage instructions are given in Appendix 
D) to carry out these calculations. 



APPENDIX B: CALCULATION OF THE 
DIFFRACTIVE TIME-SCALE AND 
BANDWIDTH 

The structure function of DM variations can be related to 
the diffractive time (jd) using Equation [5] Normally, the 
parameters of diffractive interstellar scintillation (time-scale 
Td and decorrelation frequency scale Vd) are obtained using 



DM Variations and Pulsar Timing 15 



a two dimensional auto-correlation function (ACF) of the 
dynamic spectrum S(v,t) as 

C(Au, At) = 

Np ( A \ M) A^t)A5(, + A,,t + At) (Bl) 

where N p is the number of pairs. A5*(^, t) = S(v, t) — S where 
S(v, t) is the flux density and S is the mean flux density for 
the whole observation. The diffractive parameters are de- 
fined by C(0,r d ) = C(0,0)/e and C(u d ,0) = C(0,0)/2. The 
parameters Td and Vd are obtained by fitting a 2 dimensional 
Gaussian to C(Av, At). However we often have an observed 
dynamic spectrum which is not much longer than Td and 
wider than Vd- 

We use a method based on the structure function in- 
stead of auto-correlation function to estimate the diffractive 
time-scale. In our data, the ACF is biased because we have 
few scintles in the dynamic spectra. For such cases the struc- 
ture function, defined as D(At) = [S(t) - S(t + At)] 2 /N p , 
is a better estimator because it does not require estimation 
of S. If C(At) exists, then D(At) = 2[C(0) - C(At)]. We 
estimate D(At) as 

D(At) = ^-^y ]T l AS (»> *) - AS (»> t + At )] 2 ■ ( B2 ) 

Because there is receiver noise which is white, the mea- 
sured structure function is 

D m (At) = D(At) + D W = 2C(0) - 2C{At) + D w (At) (B3) 
where D w is the structure function for the white noise. 
D w (At) = 2a 2 (At>0) 



= 0(At = 0) 



(B4) 



where a is the standard deviation of S(v,t) measured over 
the entire data-span. 

If we normalise the flux density, then S = 1. For our 
observations, the diffractive scintillation is strong and the 
observing time is much less than the refractive time-scale, 
so C(0) = 1. From the Kolmogorov spectrum, C(At) = 
exp(-(At/r d ) 5 / 3 ). 

So finally we can write the measured structure function 



D m (At) = 2[1 - exp(-(At/rd 



,5/3s 



+ D w (At) 



(B5) 



The uncertainty ao m (At) is estimated as ao m (At) = 
D m (At)y / At/T , where T is the observation time. We 
choose the equal log time interval points to fit because when 
At is large, the points are not independent. Then we can fit 
the parameters Td and D w in Equation IB5I to obtain the 
diffractive time-scale r d . 

A similar analysis can be used to obtain the diffractive 
bandwidth, v d - However, in contrast to the determination 
of Td, we do not know the form of C(Av) and, hence, it is 
not possible to fit for Vd and D m (Av). However, D w (Av) = 
2a 2 ~ D m {Avm), where Av m is the minimum frequency lag 
(in our data, it is typically 0.5 MHz). D w (Av) is the bias 
term which must be subtracted. This leads to the structure 
function being, 

D s (Av) = D m (Au) - D m {Av m ) = 2C(0) - 2C{Au). (B6) 



After normalisation (C(0) = 
tion of v d (C{v d ) = C(0)/2 
D s (Av) = 1. 



1) and according to the defini- 
= 1/2), we can obtain Vd when 



APPENDIX C: CALCULATING THE OPTIMAL 
SMOOTHING TIME 

Let t g i and t g i be idealised TO As that are affected by nei- 
ther noise nor DM variations. These TOAs correspond to 
frequencies v\ and V2 respectively where v\ > V2. Similarly, 
t g \ and t g 2o are observed TOAs at these frequencies which 
have been affected by noise and DM variations. So the ob- 
served TOAs are given by 



tqio — tq 



li(t) 



ADM(t) _ 2 



(CI) 



where rn(t) is the noise at Vi which is assumed to be white 
(i = 1,2 for the two observations respectively). The mea- 
sured estimate of DM(t) is therefore given by 

DM(i) = [t a 2o(t)-tglo(t)] _ 2 



— (tg2—tgl)— 2 



K 



ADM(t) 



-[n 2 (t) -ni(i)]- 
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(C2) 



After correcting for the DM variations, by subtracting the 
corresponding time offsets from t g \ , we obtain a set of cor- 
rected TOAs, t g ic(t), which are given by, 



tglc(t) = tglo(t) 

— tgiai — 



DM(i) _ 2 
t g 2Ci2 + ni(t)ai 



U2(t)a2 



(C3) 



and «2 = 



2 // —2 — 2\ 



where a\ — v~ /{y~ — ( 
Note that, since a\ — a,2 = 1, contributions to the tim- 
ing residuals which are frequency independent and thus the 
same in t g i and t g 2, appear unchanged in t g i c . This is an 
important property of the correction algorithm because in- 
teresting contributions such as the signature of gravitational 
waves, planets orbiting the pulsar, or ephemeris errors are 
unaltered by the correction. 

Comparing Equations IC1I and IC3I we see that the vari- 
ance of the white noise has increased from a% 1 to a% 1 a 2 + 
a N2 a 2 although the DM noise has been eliminated. We can 
improve the variance in t g i c by smoothing DM(t) before sub- 
tracting it from tgi , because smoothing reduces the white 
noise more than the DM(t) variations. The white variance 
in DM(t) is reduced by the smoothing number N s . It is 
hard to calculate the effect of smoothing DM(t) analytically, 
but we found numerically that the variance of equally sam- 
pled [DM(i)-smoothed(DM(t))] w 0.5L> D M(T sm /27r) where 
T sm — (N s — 1)to and to is the sampling time (typically 
15 d). However, the white noise in the smoothed DM(t) is 
partially correlated with that in t g i , so we have finally the 
variance in t g i c 
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We minimise this numerically. In fact, the data are not 
equally sampled so we first interpolate the raw data on to 
an equally spaced array before smoothing. To check the ef- 
fect of re-sampling we compared the theory above with a 
simulation. We realised 50 samples of DM(i) from a popula- 
tion matching DM(r) with the actual data sampling. Then 
we interpolated the simulated data onto an equally spaced 
array and found the N s which best corrected for DM(t). In 
all cases the minimum is very broad so it makes little dif- 
ference whether one uses the theoretical or simulated value 
of N s . However, the simulation is easy to implement and 
will be correct even in the case of an unusual distribution of 
samples. 



APPENDIX D: AVAILABLE SOFTWARE 

The tempo2 software was designed to allow for easy 
addition of new features and functionality in the form 
of plug-ins to the main package. During this work 
we have produced the following new plug-in packages 
which are now available as part of the tempo2 dis- 
tribution (full details are available from our web-site, 
|http: //www, atnf . csiro . au/research/pulsar/tempo2 ): 

• CALCDM: this plug-in contains the algorithms de- 
scribed in this paper. This plug-in allows the user to calcu- 
late and plot ADM(i) and obtain the corresponding struc- 
ture function. 

• SF: calculates and plots the structure function of the 
timing residuals. 

• SImISM: allows data-sets to be simulated in order to 
study the effect of a Kolmogorov process on pulsar timing 
residuals. 

The program diffTime can be used to calculate Td 
and Vd for most pulsar observations. This software is 
available as part of the PSRchive software distribution 
(http : //psr chive . sourcef orge .net ). 

Software to simulate the effect of refraction on Td 
estimates is in the SIM 2.0 distribution from UCSD: 
http : //typhoon. ucsd. edu/~ coles/ sim2 . 0/sim2 . .html 



